Data Analysis - Conflict Project
David Leydet
2022-10-22
This project is an exploration of climate change and conflict in Pakistan. This project uses the following datasets:
Pakistan Conflict Data - Acquired from the Armed Conflict Location and Event Data Project (ACLED) on 21SEP2022. See https://acleddata.com/terms-of-use/
Regional Climate Data - ERA5 Dataset with monthly averaged temperature and precipitation data for January and July beginning in 2010 and ending in 2022. The spatial resolution is approximately 9 kilometers. Temperature data is calculated at 2m above the land surface. Precipitation data is total precipitation between forecast steps (monthly). Reanalysis model. See https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land-monthly-means?tab=overview
Log Gross National Income Per Capital Variable - From: https://globaldatalab.org/shdi/metadata/lgnic/
Population Data acquired on 07NOV2022 from: https://datacommons.org/place/country/PAK?utm_medium=explore&mprop=count&popt=Person&hl=en
Pakistan Gross Domestic Product acquired on 07NOV2022 from: https://data.worldbank.org/indicator/NY.GDP.MKTP.CD?end=2021&locations=PK&start=2010
Initial Data Setup
## Set working directory
setwd("~/Desktop/University of Utah PhD /Research/r_code")
## Load packages we are using
library(mapview) #mapping package
library(raster) #raster data manipulation (Climate Data)
library(RColorBrewer) #color palettes for visualization
library(sf) #simple features for spatial data
library(tmap) #mapping package
library(viridis) #color palette for visualization
library(ncdf4) #working with netCDF files (Climate Data)
library(leaflet) #basemaps for mapview
library(ggplot2) #better figures
library(ggcorrplot) #Load the correlation plot package
library(plotly) #interactive figures
library(maps) #mapping
library(kableExtra) #creating better tables and outputs
library(dplyr) #count and data functions
library(reshape2) ## Package used to reformat data - wide to long
library(tidyverse) ##Formatting dataframes, merge, and joinPakistan Conflict Data
## Initial
## Pakistan Conflict Data
pak.con = read.csv("../data/conflict/pak_conflict.csv")
## Change data to a factor
#year don't set as a factor
## pak.con$year = factor(pak.con$year)
#interaction - this is the interaction code between actors (see ACLED codebook)
pak.con$interaction = factor(pak.con$interaction)
#inter1 - this is the type of actor for actor 1 (see ACLED codebook)
pak.con$inter1 = factor(pak.con$inter1)
#inter2 - this is the type of actor for actor 2(see ACLED codebook)
pak.con$inter2 = factor(pak.con$inter2)## Read in country shape file
countries = st_read("../data/conflict/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp")## Reading layer `ne_50m_admin_0_countries' from data source
## `/Users/davidleydet/Desktop/University of Utah PhD /Research/data/conflict/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 241 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS: WGS 84## Read in administrative boundary shape file
pak.bound = st_read("../data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm1_wfp_20220909.shp")## Reading layer `pak_admbnda_adm1_wfp_20220909' from data source
## `/Users/davidleydet/Desktop/University of Utah PhD /Research/data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm1_wfp_20220909.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 60.8786 ymin: 23.69468 xmax: 77.83397 ymax: 37.08942
## Geodetic CRS: WGS 84pak.bound.adm2 = st_read("../data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm2_wfp_20220909.shp")## Reading layer `pak_admbnda_adm2_wfp_20220909' from data source
## `/Users/davidleydet/Desktop/University of Utah PhD /Research/data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm2_wfp_20220909.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 160 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 60.8786 ymin: 23.69468 xmax: 77.83397 ymax: 37.08942
## Geodetic CRS: WGS 84## Set as a spatial feature for visualization
pak.con.sf = st_as_sf(pak.con,
coords = c("longitude", "latitude"),
crs = 4326) #WGS84## Using grepl to subset based on keywords
## keywords: water, rain, snow, river, flood
## ArcPro's output for the definition query is ~4400. ***What is the delta??***
pak.con.water.filter = dplyr::filter(pak.con.sf, grepl("water|flood|rain|river|snow", pak.con$notes, ignore.case = TRUE))
str(pak.con.water.filter)## Classes 'sf' and 'data.frame': 4761 obs. of 22 variables:
## $ data_id : int 9491475 9491460 9491461 9491470 9491741 9491742 9491743 9491775 9491394 9491739 ...
## $ event_id_cnty : chr "PAK77601" "PAK77593" "PAK77558" "PAK77585" ...
## $ event_date : chr "15-Sep-22" "13-Sep-22" "13-Sep-22" "13-Sep-22" ...
## $ year : int 2022 2022 2022 2022 2022 2022 2022 2022 2022 2022 ...
## $ event_type : chr "Protests" "Protests" "Protests" "Protests" ...
## $ sub_event_type: chr "Peaceful protest" "Peaceful protest" "Peaceful protest" "Peaceful protest" ...
## $ actor1 : chr "Protesters (Pakistan)" "Protesters (Pakistan)" "Protesters (Pakistan)" "Protesters (Pakistan)" ...
## $ assoc_actor_1 : chr "" "" "" "" ...
## $ inter1 : Factor w/ 8 levels "1","2","3","4",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ actor2 : chr "" "" "" "" ...
## $ assoc_actor_2 : chr "" "" "" "" ...
## $ inter2 : Factor w/ 9 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ interaction : Factor w/ 43 levels "10","11","12",..: 36 36 36 36 36 36 36 36 36 36 ...
## $ admin1 : chr "Sindh" "Balochistan" "Sindh" "Sindh" ...
## $ admin2 : chr "Badin" "Sohbatpur" "Sujawal" "Umerkot" ...
## $ admin3 : chr "Badin" "Sohbatpur" "Mirpur Bathoro" "Pithoro" ...
## $ location : chr "Talhar" "Sohbatpur" "Mirpur Bathoro" "Shadi Palli" ...
## $ source : chr "Daily Regional Times (Pakistan)" "Daily Regional Times (Pakistan)" "Daily Regional Times (Pakistan)" "Daily Regional Times (Pakistan)" ...
## $ source_scale : chr "Subnational" "National" "Subnational" "Subnational" ...
## $ notes : chr "On 15 September 2022, residents held a protest demonstration in Talhar town (Badin, Sindh), demanding drainage "| __truncated__ "On 13 September 2022, flood affectees held a protest demonstration in Sohbatpur town (Sohbatpur, Sindh), demand"| __truncated__ "On 13 September 2022, flood affectees held a protest demonstration in Mirpur Bathoro town (Sujawal, Sindh), dem"| __truncated__ "On 13 September 2022, flood affectees held a protest demonstration in Shadi Palli town (Umerkot, Sindh), demand"| __truncated__ ...
## $ fatalities : int 0 0 0 0 0 0 0 0 0 0 ...
## $ geometry :sfc_POINT of length 4761; first list element: 'XY' num 68.8 24.9
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:21] "data_id" "event_id_cnty" "event_date" "year" ...pak.con.water.filter.sf = st_as_sf(pak.con.water.filter)
st_write(pak.con.water.filter.sf,
dsn = "./figures/pak_water_con_sf.shp",
delete_layer = TRUE) ##Overrite the existing fileInitial Visualization Exploration
## Total events
ggplot(pak.con.sf) +
geom_sf(color = "firebrick3", alpha = 0.4)
con.total.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.sf) +
tm_symbols(col = "event_type",
alpha = 0.4,
size = 0.2,
title.col = "Conflicts By Event Type") +
tm_layout(main.title = "Pakistan Conflicts: 2010 - 2022",
legend.outside = TRUE)
con.total.map
pdf("./figures/pak_conflict_total")
con.total.map
dev.off()con.water.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.water.filter) +
tm_symbols(col = "event_type",
alpha = 0.4,
size = 0.2,
title.col = "Water Conflicts By Event Type") +
tm_layout(main.title = "Pakistan Water Conflicts: 2010 - 2022",
legend.outside = TRUE)
con.water.map
water.con.by.year.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.water.filter) +
tm_facets(by = "year") +
tm_symbols(col = "dodgerblue1",
alpha = 0.5,
size = 0.2) +
tm_layout(main.title = "Pakistan Water Conflicts: 2010 - 2022",
main.title.size = 0.75,
legend.outside = FALSE,
legend.position = c("right", "bottom"))
water.con.by.year.map
pdf("./figures/pak_water_con_by_year")
water.con.by.year.map
dev.off()con.by.year.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.sf) +
tm_facets(by = "year") +
tm_symbols(col = "firebrick3",
alpha = 0.5,
size = 0.2) +
tm_layout(main.title = "Pakistan Conflicts: 2010 - 2022",
legend.outside = TRUE)
con.by.year.map
pdf("./figures/pak_con_by_year")
con.by.year.map
dev.off()# Subset to exclude 0 fatality events
pak.con.sf.fat = subset(pak.con.sf, fatalities > 10)
## NOT WORKING PROPERLY*****
tm_shape(countries , bbox = st_bbox(pak.con.sf.fat)) +
tm_borders(col = "gray") +
tm_shape(pak.con.sf.fat) +
tm_symbols(col = "fatalities",
palette = "Reds",
alpha = 0.4,
title.col = "Number of Fatalities",
size = 0.2) +
tm_layout(main.title = "Pakistan Conflicts: 2010 - 2022",
legend.outside = TRUE) 
tm_shape(pak.bound) +
tm_borders() +
tm_fill(col = "ADM1_EN",
title = "Pakistan Provinces")
## Generate a new color palette
mycol = turbo(n = 8, begin = 0.5, end = 1, direction = 1)
con.by.district =
tm_shape(pak.bound) +
tm_borders() +
tm_fill(col = "ADM1_EN",
title = "Pakistan Provinces") +
tm_shape(pak.con.sf.fat) +
tm_symbols(col = "fatalities",
palette = mycol,
alpha = 0.4,
title.col = "Number of Fatalities",
size = 0.2) +
tm_layout(main.title = "Pakistan Conflicts by Province: 2010 - 2022",
legend.outside = TRUE) +
tm_compass(position = c("left", "top"))
con.by.district
pdf("./figures/con_by_district")
con.by.district
dev.off()tm_shape(pak.bound.adm2) +
tm_borders() +
tm_fill(col = "ADM2_EN",
title = "Pakistan Divisions") +
tm_layout(legend.show = FALSE)
Initial Data Analysis
## Create the counts by year
year.2022 = as.numeric(nrow(filter(pak.con, year == 2022)))
year.2021 = as.numeric(nrow(filter(pak.con, year == 2021)))
year.2020 = as.numeric(nrow(filter(pak.con, year == 2020)))
year.2019 = as.numeric(nrow(filter(pak.con, year == 2019)))
year.2018 = as.numeric(nrow(filter(pak.con, year == 2018)))
year.2017 = as.numeric(nrow(filter(pak.con, year == 2017)))
year.2016 = as.numeric(nrow(filter(pak.con, year == 2016)))
year.2015 = as.numeric(nrow(filter(pak.con, year == 2015)))
year.2014 = as.numeric(nrow(filter(pak.con, year == 2014)))
year.2013 = as.numeric(nrow(filter(pak.con, year == 2013)))
year.2012 = as.numeric(nrow(filter(pak.con, year == 2012)))
year.2011 = as.numeric(nrow(filter(pak.con, year == 2011)))
year.2010 = as.numeric(nrow(filter(pak.con, year == 2010)))
## Create a new data frame for the counts
# vector of years
years = c("2010", "2011", "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019", "2020", "2021", "2022")
# vector of years as a number for plotting purposes
years.num = seq(from = 2010, to = 2022, by = 1)
# vector of conflict counts
num.conflicts = c(year.2010, year.2011, year.2012, year.2013, year.2014, year.2015, year.2016, year.2017, year.2018, year.2019, year.2020, year.2021, year.2022)
# Combine the vectors into a data frame
con.by.year = data.frame (years, num.conflicts)
# Combine the vectors using the years as a number
con.by.year.num = data.frame(years.num, num.conflicts)
# Visualize the Count
# How does one add a trendline?
num.conflicts.sctrplot =
ggplot(data = con.by.year.num, aes(x = years.num, y = num.conflicts)) +
geom_point(shape = 18, size = 3, color = "red") +
geom_smooth(method = lm,
se = FALSE) +
xlab("Year") +
ylab("Number of Conflicts") +
labs(title = "Pakistan Conflicts by Year") +
scale_x_discrete(limits = c(years.num)) +
theme_bw()
num.conflicts.sctrplot## `geom_smooth()` using formula 'y ~ x'
pdf("./figures/con_by_year_scatterplot")
num.conflicts.sctrplot
dev.off()##Count the number of water conflicts by year
water.year.2022 = as.numeric(nrow(filter(pak.con.water.filter, year == 2022)))
water.year.2021 = as.numeric(nrow(filter(pak.con.water.filter, year == 2021)))
water.year.2020 = as.numeric(nrow(filter(pak.con.water.filter, year == 2020)))
water.year.2019 = as.numeric(nrow(filter(pak.con.water.filter, year == 2019)))
water.year.2018 = as.numeric(nrow(filter(pak.con.water.filter, year == 2018)))
water.year.2017 = as.numeric(nrow(filter(pak.con.water.filter, year == 2017)))
water.year.2016 = as.numeric(nrow(filter(pak.con.water.filter, year == 2016)))
water.year.2015 = as.numeric(nrow(filter(pak.con.water.filter, year == 2015)))
water.year.2014 = as.numeric(nrow(filter(pak.con.water.filter, year == 2014)))
water.year.2013 = as.numeric(nrow(filter(pak.con.water.filter, year == 2013)))
water.year.2012 = as.numeric(nrow(filter(pak.con.water.filter, year == 2012)))
water.year.2011 = as.numeric(nrow(filter(pak.con.water.filter, year == 2011)))
water.year.2010 = as.numeric(nrow(filter(pak.con.water.filter, year == 2010)))
##Automate this: function script?
##Function Works!
##Automate loop??
mycountbyyear = function(x) {
as.numeric(nrow(filter(x, year == 2010)))
}
testyear = mycountbyyear(pak.con.water.filter)
##Create a vector of water conflicts by year
water.num.conflicts = c(water.year.2010, water.year.2011, water.year.2012, water.year.2013, water.year.2014, water.year.2015, water.year.2016, water.year.2017, water.year.2018, water.year.2019, water.year.2020, water.year.2021, water.year.2022)
# Combine the vectors into a data frame. Need to use years.num in order to plot a trend line
water.con.by.year = data.frame (years.num, water.num.conflicts)
##Plot it!
num.water.conflicts.sctrplot =
ggplot(data = water.con.by.year, aes(x = years.num, y = water.num.conflicts)) +
geom_point(shape = 18, size = 3, color = "blue") +
geom_smooth(method = NULL,
se = TRUE,
color = "black",
linetype = "dashed",
size = 0.5) +
xlab("Year") +
ylab("Number of Water Conflicts") +
labs(title = "Pakistan Water Conflicts by Year") +
scale_x_discrete(limits = c(years.num)) +
theme_bw() ## Warning: Continuous limits supplied to discrete scale.
## Did you mean `limits = factor(...)` or `scale_*_continuous()`?num.water.conflicts.sctrplot## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
pdf("./figures/water_con_by_year_scatterplot")
num.water.conflicts.sctrplot
dev.off()Pakistan GNI Data
## Initial Read
## ****Note - the years are from 1990 - 2019****
## Log Gross National Income per capita in thousands of US Dollars (2011 PPP)
## Ask Simon about time series data
pak.gni = read.csv("../data/conflict/pak_subnat_gdp.csv")
## Convert it from wide format to long format
pak.gni.long = melt(pak.gni, id.vars = "Region")
## Change the column name to year
colnames(pak.gni.long)[2] = "year"
## Remove the "x" from the years
pak.gni.long$year = gsub("X", "", as.factor(pak.gni.long$year))
## Subset the temporal timeframe: 2010 - 2022 (2019 in this case)
###pak.gni.time.subset = subset(pak.gni.long, Factor == "2010")pak.gni.plot =
ggplot(data = pak.gni.long, aes(x = year, y = value, color = Region)) +
geom_point() +
xlab("Year") +
ylab("GNI per Capita (Thousands USD)") +
# xlim(2010, 2019) +
theme(axis.text.x = element_text(angle = 90))
theme_bw()## List of 93
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 11
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ title : NULL
## $ aspect.ratio : NULL
## $ axis.title : NULL
## $ axis.title.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.75points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.75points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.bottom : NULL
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.75points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y.left : NULL
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.75points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.2points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.2points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.bottom : NULL
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 1
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.2points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y.left : NULL
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.2points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.x : NULL
## $ axis.ticks.x.top : NULL
## $ axis.ticks.x.bottom : NULL
## $ axis.ticks.y : NULL
## $ axis.ticks.y.left : NULL
## $ axis.ticks.y.right : NULL
## $ axis.ticks.length : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom: NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.x.top : NULL
## $ axis.line.x.bottom : NULL
## $ axis.line.y : NULL
## $ axis.line.y.left : NULL
## $ axis.line.y.right : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ legend.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.key.size : 'simpleUnit' num 1.2lines
## ..- attr(*, "unit")= int 3
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.align : NULL
## $ legend.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.title.align : NULL
## $ legend.position : chr "right"
## $ legend.direction : NULL
## $ legend.justification : chr "center"
## $ legend.box : NULL
## $ legend.box.just : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "unit")= int 1
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.box.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ panel.background :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.border :List of 5
## ..$ fill : logi NA
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.spacing : 'simpleUnit' num 5.5points
## ..- attr(*, "unit")= int 8
## $ panel.spacing.x : NULL
## $ panel.spacing.y : NULL
## $ panel.grid :List of 6
## ..$ colour : chr "grey92"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major : NULL
## $ panel.grid.minor :List of 6
## ..$ colour : NULL
## ..$ size : 'rel' num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major.x : NULL
## $ panel.grid.major.y : NULL
## $ panel.grid.minor.x : NULL
## $ panel.grid.minor.y : NULL
## $ panel.ontop : logi FALSE
## $ plot.background :List of 5
## ..$ fill : NULL
## ..$ colour : chr "white"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ plot.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.title.position : chr "panel"
## $ plot.subtitle :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : num 1
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 5.5points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption.position : chr "panel"
## $ plot.tag :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag.position : chr "topleft"
## $ plot.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ strip.background :List of 5
## ..$ fill : chr "grey85"
## ..$ colour : chr "grey20"
## ..$ size : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ strip.background.x : NULL
## $ strip.background.y : NULL
## $ strip.placement : chr "inside"
## $ strip.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey10"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 4.4points 4.4points 4.4points 4.4points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.x : NULL
## $ strip.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.switch.pad.grid : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ strip.switch.pad.wrap : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ strip.text.y.left :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUEpak.gni.plot
## Save as a .pdf
ggsave("./figures/gni_by_year_district.pdf", pak.gni.plot)##2022 GDP has not been published yet, so we estimated based on an average gdp growth of 3.7%
pak.gdp = read.csv("../data/conflict/pak_gdp_2010_2022.csv")
pak.gdp.plot =
ggplot(data = pak.gdp, aes(x = year, y = gdp)) +
geom_point(col = "forestgreen") +
geom_line(linetype = "dashed",
col = "forestgreen") +
labs(title = "Pakistan GDP (2010-2022)") +
xlab("Year") +
ylab("GDP USD($)")
pak.gdp.plot
#May want to scale the GDPpdf("./figures/pak_gdp_by_year")
pak.gdp.plot
dev.off()Pakistan Climate Data
Temperature Data
## Read in the ERA5 climate data
pak.temp = raster("../data/conflict/pak_clim.nc", varname = "t2m")
## Look at the data
##26 bands (layers)
##13 years worth of data * 2 months(Jan, Jul) = 26 observations (bands?)
## Syntax variable = raster("filepath", band = #) to extract the specific band you need
## library(ncdf4) this package allows you to read into the metadata
## Could use stack() to look at all of the bands
pak.temp## class : RasterLayer
## band : 1 (of 26 bands)
## dimensions : 201, 251, 50451 (nrow, ncol, ncell)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : pak_clim.nc
## names : X2.metre.temperature
## z-value : 2010-01-01
## zvar : t2m## CRS Check
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs ## Temperature is in Kelvin (how do I convert the scale to Celsius? Kelvin Temp - 273.15)
## January Temperature. How to plot July?
mytemp.pal = brewer.pal(n = 9, name = "OrRd")
plot(pak.temp,
main = "January Temperature - 2010",
col = mytemp.pal)
plot(st_geometry(pak.bound),
border = "black",
add = TRUE)
##Not working correctly now? 20221107
pak.temp.masked = mask(pak.temp, mask = pak.bound)
plot(pak.temp.masked,
main = "Pakistan Jan 2010 Temperature - Masked",
col = mytemp.pal)
plot(st_geometry(pak.bound), add = TRUE)
##Use stack() to load all of the temperature bands
temp.stack = stack("../data/conflict/pak_clim.nc", varname = "t2m")
##View the stack data
temp.stack## class : RasterStack
## dimensions : 201, 251, 50451, 26 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2010.01.01, X2010.07.01, X2011.01.01, X2011.07.01, X2012.01.01, X2012.07.01, X2013.01.01, X2013.07.01, X2014.01.01, X2014.07.01, X2015.01.01, X2015.07.01, X2016.01.01, X2016.07.01, X2017.01.01, ...
## Date/time : 2010-01-01 - 2022-07-01 (range)##January 2010 stack
temp.stack.jan2010 = temp.stack[[1]]
temp.stack.jan2010## class : RasterLayer
## band : 1 (of 26 bands)
## dimensions : 201, 251, 50451 (nrow, ncol, ncell)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : pak_clim.nc
## names : X2010.01.01
## Date/time : 2010-01-01
## zvar : t2m##Crop the temperature stack by the Pakistan Boundary
temp.stack.crop = crop(temp.stack, pak.bound)
temp.stack.crop## class : RasterBrick
## dimensions : 134, 170, 22780, 26 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1 (x, y)
## extent : 60.85, 77.85, 23.65, 37.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : X2010.01.01, X2010.07.01, X2011.01.01, X2011.07.01, X2012.01.01, X2012.07.01, X2013.01.01, X2013.07.01, X2014.01.01, X2014.07.01, X2015.01.01, X2015.07.01, X2016.01.01, X2016.07.01, X2017.01.01, ...
## min values : 242.9010, 262.5190, 239.7853, 264.1235, 238.9400, 261.9992, 239.1240, 262.7146, 240.7158, 262.4167, 240.7767, 265.0068, 240.8986, 265.3265, 242.3812, ...
## max values : 294.2676, 310.3270, 293.8248, 310.1821, 293.5142, 311.2782, 293.4889, 311.3368, 292.2054, 311.5416, 293.4717, 311.2161, 295.2015, 311.4714, 293.9168, ...
## time : 2010-01-01, 2022-07-01 (min, max)##Color Template for Temperature
##Create a function to add the Pakistan Level 2 boundaries to the figure
addborder = function(){
plot(as_Spatial(pak.bound),
add = TRUE)
}
plot(temp.stack.crop[[1]],
col = mytemp.pal,
zlim = c(250, 320), ##zlim makes sure the scale is the same
addfun = addborder)
pak.temp.stack.mask = mask(temp.stack, mask = pak.bound)
plot(pak.temp.stack.mask,
col = mytemp.pal,
zlim = c(250, 320), ##zlim makes sure the scale is the same
addfun = addborder)
##Use the temperature stack (masked) to calculate the mean average temperature for the country
temp.avg.all = cellStats(pak.temp.stack.mask, mean)
##Visualize the data - this differs slighlt from running the entire raster data (without the mask)
#plot(temp.avg.all,
#type = "p")
temp.avg.all## X2010.01.01 X2010.07.01 X2011.01.01 X2011.07.01 X2012.01.01 X2012.07.01
## 283.0038 300.5736 281.0138 300.7901 279.7675 301.7816
## X2013.01.01 X2013.07.01 X2014.01.01 X2014.07.01 X2015.01.01 X2015.07.01
## 280.9438 301.1685 280.7888 301.3565 281.6683 300.4983
## X2016.01.01 X2016.07.01 X2017.01.01 X2017.07.01 X2018.01.01 X2018.07.01
## 282.7455 301.3774 281.1049 300.9322 282.4930 301.4175
## X2019.01.01 X2019.07.01 X2020.01.01 X2020.07.01 X2021.01.01 X2021.07.01
## 281.5669 301.5670 279.5454 301.7934 280.6591 301.7339
## X2022.01.01 X2022.07.01
## 280.8029 299.3014##Shape the average temperature values into a data frame
temp.avg.df = as.data.frame(temp.avg.all)
temp.avg.df## temp.avg.all
## X2010.01.01 283.0038
## X2010.07.01 300.5736
## X2011.01.01 281.0138
## X2011.07.01 300.7901
## X2012.01.01 279.7675
## X2012.07.01 301.7816
## X2013.01.01 280.9438
## X2013.07.01 301.1685
## X2014.01.01 280.7888
## X2014.07.01 301.3565
## X2015.01.01 281.6683
## X2015.07.01 300.4983
## X2016.01.01 282.7455
## X2016.07.01 301.3774
## X2017.01.01 281.1049
## X2017.07.01 300.9322
## X2018.01.01 282.4930
## X2018.07.01 301.4175
## X2019.01.01 281.5669
## X2019.07.01 301.5670
## X2020.01.01 279.5454
## X2020.07.01 301.7934
## X2021.01.01 280.6591
## X2021.07.01 301.7339
## X2022.01.01 280.8029
## X2022.07.01 299.3014##Write to csv to transform....need to find a better way to do this!!!!***
write.csv(temp.avg.df,
"../data/conflict/pak_avg_temp")##Read in the converted csv
pak.temp.avg = read.csv("../data/conflict/pak_avg_temp_converted.csv")
pak.temp.avg## year avg_jan_temp avg_jul_temp avg_jan_celsius avg_jul_celsius
## 1 2010 283.0038 300.5736 9.853772 27.42361
## 2 2011 281.0138 300.7901 7.863778 27.64015
## 3 2012 279.7675 301.7816 6.617535 28.63159
## 4 2013 280.9438 301.1685 7.793822 28.01848
## 5 2014 280.7888 301.3565 7.638813 28.20651
## 6 2015 281.6683 300.4983 8.518308 27.34832
## 7 2016 282.7455 301.3774 9.595462 28.22739
## 8 2017 281.1049 300.9322 7.954915 27.78219
## 9 2018 282.4930 301.4175 9.342984 28.26746
## 10 2019 281.5669 301.5670 8.416871 28.41698
## 11 2020 279.5454 301.7934 6.395355 28.64344
## 12 2021 280.6591 301.7339 7.509081 28.58392
## 13 2022 280.8029 299.3014 7.652921 26.15137avg.temp.plot =
ggplot(data = pak.temp.avg) +
geom_line(aes(x=year, y = avg_jan_celsius), col = "blue") +
geom_line(aes(x=year, y = avg_jul_celsius), col = "red") +
xlab("Year") +
ylab("Temperature (Celsius)")
avg.temp.plot
Precipitation Data
## Read in the ERA5 climate data
## Value is in meters
##26 bands (layers)
##13 years worth of data * 2 months(Jan and Jul) = 26 observations (bands?)
## Syntax variable = raster("filepath", band = #) to extract the specific band you need
## library(ncdf4) this package allows you to read into the metadata
## Could use stack() to look at all of the bands
pak.precip = raster("../data/conflict/pak_clim.nc", varname = "tp")
## Look at the data
pak.precip## class : RasterLayer
## band : 1 (of 26 bands)
## dimensions : 201, 251, 50451 (nrow, ncol, ncell)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : pak_clim.nc
## names : Total.precipitation
## z-value : 2010-01-01
## zvar : tp##January precipitation
myprecip.pal = brewer.pal(n = 9, name = "Blues")
plot(pak.precip,
main = "January Precipitation - 2010",
col = myprecip.pal)
plot(st_geometry(pak.bound.adm2),
border = "black",
add = TRUE)
##Use stack() to load all of the precipitation bands
##REMINDER This is total precipitation in meters
precip.stack = stack("../data/conflict/pak_clim.nc", varname = "tp")
##View the stack data
precip.stack## class : RasterStack
## dimensions : 201, 251, 50451, 26 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2010.01.01, X2010.07.01, X2011.01.01, X2011.07.01, X2012.01.01, X2012.07.01, X2013.01.01, X2013.07.01, X2014.01.01, X2014.07.01, X2015.01.01, X2015.07.01, X2016.01.01, X2016.07.01, X2017.01.01, ...
## Date/time : 2010-01-01 - 2022-07-01 (range)pak.precip.stack.mask = mask(precip.stack, mask = pak.bound)
plot(pak.precip.stack.mask,
col = myprecip.pal,
zlim = c(0, 0.0030), ##zlim makes sure the scale is the same
addfun = addborder)
##Use the precipitation stack (masked) to calculate the mean average temperature for the country
precip.avg.all = cellStats(pak.precip.stack.mask, mean)
##Visualize the data - this differs slighlt from running the entire raster data (without the mask)
#plot(precip.avg.all,
#type = "p")
precip.avg.all## X2010.01.01 X2010.07.01 X2011.01.01 X2011.07.01 X2012.01.01 X2012.07.01
## 0.0005839226 0.0043177244 0.0003963835 0.0020611769 0.0009019366 0.0010772682
## X2013.01.01 X2013.07.01 X2014.01.01 X2014.07.01 X2015.01.01 X2015.07.01
## 0.0003545558 0.0018138459 0.0003568661 0.0015595845 0.0007123514 0.0038597275
## X2016.01.01 X2016.07.01 X2017.01.01 X2017.07.01 X2018.01.01 X2018.07.01
## 0.0006265805 0.0022160874 0.0019144569 0.0026588202 0.0002786720 0.0019427409
## X2019.01.01 X2019.07.01 X2020.01.01 X2020.07.01 X2021.01.01 X2021.07.01
## 0.0012742865 0.0020474729 0.0018617504 0.0016054658 0.0003895901 0.0021279386
## X2022.01.01 X2022.07.01
## 0.0017411640 0.0063039410##Shape the average temperature values into a data frame
precip.avg.df = as.data.frame(precip.avg.all)
precip.avg.df## precip.avg.all
## X2010.01.01 0.0005839226
## X2010.07.01 0.0043177244
## X2011.01.01 0.0003963835
## X2011.07.01 0.0020611769
## X2012.01.01 0.0009019366
## X2012.07.01 0.0010772682
## X2013.01.01 0.0003545558
## X2013.07.01 0.0018138459
## X2014.01.01 0.0003568661
## X2014.07.01 0.0015595845
## X2015.01.01 0.0007123514
## X2015.07.01 0.0038597275
## X2016.01.01 0.0006265805
## X2016.07.01 0.0022160874
## X2017.01.01 0.0019144569
## X2017.07.01 0.0026588202
## X2018.01.01 0.0002786720
## X2018.07.01 0.0019427409
## X2019.01.01 0.0012742865
## X2019.07.01 0.0020474729
## X2020.01.01 0.0018617504
## X2020.07.01 0.0016054658
## X2021.01.01 0.0003895901
## X2021.07.01 0.0021279386
## X2022.01.01 0.0017411640
## X2022.07.01 0.0063039410##Write to csv to transform....need to find a better way to do this!!!!***
write.csv(precip.avg.df,
"../data/conflict/pak_avg_precip")##Read in the converted csv
pak.precip.avg = read.csv("../data/conflict/pak_avg_precip_converted.csv")
pak.precip.avg## year avg_jan_precip avg_jul_precip avg_jan_precip_mm avg_jul_precip_mm
## 1 2010 0.000583923 0.004317724 0.5839226 4.317724
## 2 2011 0.000396383 0.002061177 0.3963835 2.061177
## 3 2012 0.000901937 0.001077268 0.9019366 1.077268
## 4 2013 0.000354556 0.001813846 0.3545558 1.813846
## 5 2014 0.000356866 0.001559585 0.3568661 1.559585
## 6 2015 0.000712351 0.003859728 0.7123514 3.859728
## 7 2016 0.000626581 0.002216087 0.6265805 2.216087
## 8 2017 0.001914457 0.002658820 1.9144569 2.658820
## 9 2018 0.000278672 0.001942741 0.2786720 1.942741
## 10 2019 0.001274286 0.002047473 1.2742865 2.047473
## 11 2020 0.001861750 0.001605466 1.8617504 1.605466
## 12 2021 0.000389590 0.002127939 0.3895901 2.127939
## 13 2022 0.001741164 0.006303941 1.7411640 6.303941avg.precip.plot =
ggplot(data = pak.precip.avg) +
geom_line(aes(x=year, y = avg_jan_precip_mm, color = "avg_jan_precip")) +
geom_line(aes(x=year, y = avg_jul_precip_mm, color = "avg_jul_precip")) +
xlab("Year") +
ylab("Total Precipitation (mm)") +
scale_color_manual(name = "Seasonal Precipitation",
values = c("avg_jan_precip" = "blue",
"avg_jul_precip" = "red"))
avg.precip.plot
Population Data
##Data acquired from datacommons.org (world bank)
##2022 is estimated based on the 10-year average growth rate
pak.pop = read.csv("../data/conflict/pak_pop.csv")
ggplot(data = pak.pop, aes(x = year, y = pak_pop)) +
geom_point(col = "darkorange2") +
geom_line(linetype = "dashed",
col = "darkorange2") +
labs(title = "Pakistan Population (2010-2022)") +
xlab("Year") +
ylab("Population")
Models Build
##Test: Merge the GDP and Number of Conflicts together
##Change the column header in order to merge
colnames(con.by.year)[1] = "year"
##Merge
pak.model.df = merge(con.by.year, pak.gdp, by = "year")
##It worked!
##Add water conflicts
##use by.x and by.y to merge using different variable names
pak.model.df = merge(pak.model.df, water.con.by.year, by.x ="year", by.y = as.character("years.num"))
##Now to add population
pak.model.df = merge(pak.model.df, pak.pop, by = "year")
##Convert year to number
pak.model.df$year = as.numeric(pak.model.df$year)
##Merge temperature data
pak.model.df = merge(pak.model.df, pak.temp.avg, by = "year")
##Merge precipitation data
pak.model.df = merge(pak.model.df, pak.precip.avg, by = "year")##Write to csv to save***
write.csv(pak.model.df,
"../data/conflict/pak_model_df")#pairs(pak.model.df)
ggcorrplot(cor(pak.model.df),
method = "square",
type = "lower",
lab = TRUE,
colors = c("blue", "darksalmon", "firebrick"))
#p.mat = cor_pmat(pak.model.df)) #x out non-significant p-valuesGLM - Poisson
##Starting with a generalized linear model - Poisson - as the number of conflicts is count data
pak.nat.glm = glm(num.conflicts ~ year + gdp + pak_pop + avg_jan_celsius + avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm,
family = poisson(link = 'log'),
data = pak.model.df)
summary(pak.nat.glm)##
## Call:
## glm(formula = num.conflicts ~ year + gdp + pak_pop + avg_jan_celsius +
## avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm,
## family = poisson(link = "log"), data = pak.model.df)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.5569 -9.4196 6.0944 14.9581 -12.9770 4.7249 -13.0862 -5.8290
## 9 10 11 12 13
## 13.3055 6.3542 0.1537 -7.4487 1.2489
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.594e+03 1.418e+02 11.241 < 2e-16 ***
## year -8.140e-01 7.199e-02 -11.308 < 2e-16 ***
## gdp -7.393e-13 2.446e-13 -3.022 0.00251 **
## pak_pop 2.003e-07 1.696e-08 11.810 < 2e-16 ***
## avg_jan_celsius -3.162e-02 6.792e-03 -4.655 3.23e-06 ***
## avg_jul_celsius 5.432e-01 2.203e-02 24.661 < 2e-16 ***
## avg_jan_precip_mm 3.746e-02 8.052e-03 4.653 3.27e-06 ***
## avg_jul_precip_mm 1.655e-01 1.307e-02 12.661 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4909.8 on 12 degrees of freedom
## Residual deviance: 1020.4 on 5 degrees of freedom
## AIC: 1172.7
##
## Number of Fisher Scoring iterations: 4##Need to transform the coefficients ##Coefficient transformation as the output is in the log values (link = log)
##Reminder Interpretation Note:
##Log odds scale to odds scale - it becomes multiplicative - a rate of change (bigger than 1 it is a positive change, if it is less than 1 it is a negative change)
pak.nat.glm.coef = exp(coef(pak.nat.glm))
pak.nat.glm.coef## (Intercept) year gdp pak_pop
## Inf 0.4430800 1.0000000 1.0000002
## avg_jan_celsius avg_jul_celsius avg_jan_precip_mm avg_jul_precip_mm
## 0.9688734 1.7215653 1.0381726 1.1799515Residuals
##View the residuals
plot(residuals.glm(pak.nat.glm))
abline(h = 0, col = "red")
##Fitting a spline may be better for this model?plot(pak.nat.glm)
